﻿using System;
using System.Collections.Generic;
using System.Text;
using System.Runtime.InteropServices;
//////////////////////////////////////////////////////////////////////////
//Writen by Liu Yonghe,Henan Polytechnic University
//河南理工大学刘永和编写
//////////////////////////////////////////////////////////////////////////
namespace GeoFly
{
    public class CUDA_PitFill
    {
        /// <summary>
        /// 普通串行算法
        /// </summary>
        /// <param name="DirectionMatrix"></param>
        /// <param name="dem"></param>
        /// <returns></returns>
        public double[,] CheckPit_SweepLine(int[,] DirectionMatrix, GridLayer dem)
        {
            double[,] DEM = dem.Values;
            int[,]  marks = new int[dem.rowCount, dem.colCount];
            double[,] waterData = new double[dem.rowCount, dem.colCount];
            for (int row = 0; row < dem.rowCount; row++)
            {
                waterData[row,0]=DEM[row,0];
                waterData[row,dem.colCount-1]=DEM[row,dem.colCount-1];
            }
            for (int col = 0; col < dem.colCount; col++)
            {
                waterData[0, col] = DEM[0, col];
                waterData[dem.rowCount - 1, col] = DEM[dem.rowCount - 1, col];
            }
            for (int row = 1; row < dem.rowCount - 1; row++)
            {
                for (int col = 1; col < dem.colCount - 1; col++)
                    waterData[row, col] = 1e20;
            }
           
            double epsilon = 1e-10;
            //从左到右、从下至上扫描
            for (int time = 0; time <3; time++)
            {
                for (int col = 1; col < dem.colCount - 1; col++)
                {
                    for (int row = 1; row < dem.rowCount - 1; row++)
                    {
                        //获取前一行的三个邻域像元和当前列的下方像元
                        int[] r = new int[] { row - 1, row, row + 1, row - 1 };
                        int[] c = new int[] { col - 1, col - 1, col - 1, col };
                        //找出邻域中高程最低的水面

                        double min = 1e20;
                        for (int i = 0; i < 4; i++)
                        {
                            int row_n = r[i];
                            int col_n = c[i];
                            double neighbor = waterData[row_n, col_n];
                            if (min > neighbor)
                                min = neighbor;
                        }
                        double demv = DEM[row, col];
                       
                        double value = min > demv ? min : demv;// +epsilon;
                        if (value < waterData[row, col])
                            waterData[row, col] = value + epsilon;
                        //double dem_v = DEM[row, col];
                    }
                }
                //从右到左、从上至下扫描
                for (int col = dem.colCount - 2; col > 0; col--)
                {
                    for (int row = dem.rowCount - 2; row > 0; row--)
                    {
                        if (marks[row, col] > 0)
                            continue;
                        //获取前一行的三个邻域像元和当前列的下方像元
                        int[] r = new int[] { row + 1, row, row - 1, row + 1 };
                        int[] c = new int[] { col + 1, col + 1, col + 1, col };

                        double min = 1e20;
                        for (int i = 0; i < 4; i++)
                        {
                            int row_n = r[i];
                            int col_n = c[i];
                            double neighbor = waterData[r[i], c[i]];
                            if (min > neighbor)
                                min = neighbor;
                        }
                        double demv = DEM[row, col];
                        double value = min > demv ? min : demv;// +epsilon;
                        if (value < waterData[row, col])
                            waterData[row, col] = value + epsilon;
                    }
                }
            }
            //从下至上扫描
            //for (int row = 1; row < dem.rowCount - 1; row++)
            //{
            //    for (int col = 1; col < dem.colCount - 1; col++)
            //    {
            //        double value= Math.Max(waterData[row - 1, col], DEM[row, col]);
            //        if(value< waterData[row, col])
            //            waterData[row,col]=value+epsilon;
            //    }
            //}
            //for (int row = dem.rowCount - 2; row > 0; row--)
            //{
            //    for (int col = dem.colCount - 2; col > 0; col--)
            //    {
            //        double value = Math.Max(waterData[row + 1, col], DEM[row, col]);
            //        if (value < waterData[row, col])
            //            waterData[row, col] = value+epsilon;
            //    }
            //}
            dem.Values = waterData;
            return waterData;
        }
        public double[,] CheckPit_SweepLine2(int[,] DirectionMatrix, GridLayer dem)
        {
            double[,] DEM = dem.Values;
            int[,] marks = new int[dem.rowCount, dem.colCount];
            double[,] waterData = new double[dem.rowCount, dem.colCount];
            for (int row = 0; row < dem.rowCount; row++)
            {
                waterData[row, 0] = DEM[row, 0];
                waterData[row, dem.colCount - 1] = DEM[row, dem.colCount - 1];
            }
            for (int col = 0; col < dem.colCount; col++)
            {
                waterData[0, col] = DEM[0, col];
                waterData[dem.rowCount - 1, col] = DEM[dem.rowCount - 1, col];
            }
            for (int row = 1; row < dem.rowCount - 1; row++)
            {
                for (int col = 1; col < dem.colCount - 1; col++)
                    waterData[row, col] = 1e20;
            }

            double epsilon = 1e-10;
            //从左到右、从下至上扫描

            for (int time = 0; time < 1; time++)
            {
                for (int col = 1; col < dem.colCount - 1; col++)
                {
                    for (int row = 1; row < dem.rowCount - 1; row++)
                    {
                        //获取前一行的三个邻域像元和当前列的下方像元
                        int[] r = new int[] { row - 1, row, row + 1 };
                        int[] c = new int[] { col - 1, col - 1, col - 1 };
                        //找出邻域中高程最低的水面

                        double min = 1e20;
                        for (int i = 0; i < 3; i++)
                        {
                            int row_n = r[i];
                            int col_n = c[i];
                            double neighbor = waterData[row_n, col_n];
                            if (min > neighbor)
                                min = neighbor;
                        }
                        double demv = DEM[row, col];
                        double value = min > demv ? min : demv;// +epsilon;
                        if (value < waterData[row, col])
                            waterData[row, col] = value + epsilon;
                        //double dem_v = DEM[row, col];
                    }
                }
                //从右到左、从上至下扫描
                for (int col = dem.colCount - 2; col > 0; col--)
                {
                    for (int row = dem.rowCount - 2; row > 0; row--)
                    {
                        if (marks[row, col] > 0)
                            continue;
                        //获取前一行的三个邻域像元和当前列的下方像元
                        int[] r = new int[] { row + 1, row, row - 1 };
                        int[] c = new int[] { col + 1, col + 1, col + 1 };

                        double min = 1e20;
                        for (int i = 0; i < 3; i++)
                        {
                            int row_n = r[i];
                            int col_n = c[i];
                            double neighbor = waterData[r[i], c[i]];
                            if (min > neighbor)
                                min = neighbor;
                        }
                        double demv = DEM[row, col];
                        double value = min > demv ? min : demv;// +epsilon;
                        if (value < waterData[row, col])
                            waterData[row, col] = value + epsilon;
                    }
                }
                //从下至上扫描
                for (int row = 1; row < dem.rowCount - 1; row++)
                {
                    for (int col = 1; col < dem.colCount - 1; col++)
                    {
                        double value = Math.Max(waterData[row - 1, col], DEM[row, col]);
                        if (value < waterData[row, col])
                            waterData[row, col] = value + epsilon;
                    }
                }
                for (int row = dem.rowCount - 2; row > 0; row--)
                {
                    for (int col = dem.colCount - 2; col > 0; col--)
                    {
                        double value = Math.Max(waterData[row + 1, col], DEM[row, col]);
                        if (value < waterData[row, col])
                            waterData[row, col] = value + epsilon;
                    }
                }
            }
            dem.Values = waterData;
            return waterData;
        }
        public float[] DEM;
       
    }
}
